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We examine the use of synchronization as a mechanism for extracting parameter and state in- 
formation from experimental systems. We focus on important aspects of this problem that have 
received little attention previously, and we explore them using experiments and simulations with the 
chaotic Colpitis oscillator as an example system. We explore the impact of model imperfection on 
the ability to extract valid information from an experimental system. We compare two optimization 
methods: an initial value method and a constrained method. Each of these involve coupling the 
model equations to the experimental data in order to regularize the chaotic motions on the syn- 
chronization manifold. We explore both time dependent and time independent coupling. We also 
examine both optimized and fixed (or manually adjusted) coupling. For the case of an optimized 
time dependent coupling function u{t) we find a robust structure which includes sharp peaks and 
intervals where it is zero. This structure shows a strong correlation with the location in phase space 
and appears to depend on noise, imperfections of the model, and the Lyapunov direction vectors. 
For time independent coupling we find the counterintuitive result that often the optimal rms er- 
ror in fitting the model to the data initially increases with coupling strength. Comparison of this 
result with that obtained using simulated data may provide one measure of model imperfection. 
The constrained method with time dependent coupling appears to have benefits in synchronizing 
long datasets with minimal impact, while the initial value method with time independent coupling 
tends to be substantially faster, more flexible and easier to use. We also describe a new method of 
coupling which is useful for sparse experimental data sets. Our use of the Colpitis oscillator allows 
us to explore in detail the case of a system with one positive Lyapunov exponent. The methods 
we explored are easily extended to driven systems such as neurons with time dependent injected 
current. They are expected to be of value in nonchaotic systems as well. Software is available on 
request. 



I. INTRODUCTION 

Physical models of nonlinear systems formulated as 
differential equations or as discrete time maps typically 
have unknown parameters representing our lack of de- 
tailed knowledge or our conjectures about how the dy- 
namics operates in these systems. Further, forecasting 
the future of such systems requires the knowledge both 
of these parameters as well as the knowledge of the full 
state of the system at the time a forecast begins. Un- 
fortunately, in interesting complex dynamics as diverse 
as weather forecasting to prediction of networks of bio- 
logical neurons, we cannot measure the full states of the 
system, and we must rely on an accurate parameter and 
state estimation procedure that allows reliable inference 
both of the parameters and of the full states of the dy- 
namical systems from observations of some subset of the 
dynamical variables. 
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We make use of thephenomenon of chaos synchroniza- 
tion [E i, S, II i, i, 0, i, i]. This idea can be applied 
to the problem of parameter estimation by treating the 
measured data as the driver system, and implementing 
the model or response system numerically. The receiver 
system is coupled to the driver system. The fixed pa- 
rameters and the initial conditions in the model/response 
system are adjusted to find the set of parameters which 
yields the smallest synchronization error of the two sys- 
tems [IS El, m E 0, [H, Ell ■ Sufficiently strong cou- 
pling can often mitigate the effects of a positive Lyapunov 
exponent in the response system. It also tends to dra- 
matically reduce the difhculty in finding a path through 
the parameter space from the initial guess to the optimal 
values. This can often mean the difference between suc- 
cess and failure, because when the number of parameters 
is relatively large, it becomes a practical impossibility to 
fully explore that space. 

In this paper we examine as a "test system" a seem- 
ingly simple experimental system, the electronic Colpitis 
Oscillator. This system has the dual advantages of being 
easy to study, while still exhibiting difficulties in model- 
ing that are commonplace in more complex systems such 
as networks of neurons. It is often assumed that a perfect 
model for a system is available, but in reality no model 
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is perfect and often the seemingly small imperfections of 
a model can have a substantial impact on our ability to 
extract accurate information about that system from the 
available data. 

In the control theory literature the task of determin- 
ing the unobserved dynamical variables from observa- 
tions of a subset of the whole collection of state variables 
goe s under the name of finding observers for the system 
[li m, m [13, mi, m . to the extent one may consider 
the unknown parameters of the system as state variables 
constant in time, the observer problem also includes the 
estimation of these parameters. 

Chaotic oscillations of the system may mean that the 
observations we make and their equivalent dynamical 
variables in the models we develop for those nonlinear 
systems may not synchronize. This will severely impede 
our ability to estimate either unknown states or param- 
eters for the system of interest. This issue was first ad- 
dressed by So, Ott, and Dayaswana [l^ in 1994 where 
they recognized that the states of the observed system, 
represented by vectors x(t), and the states of the model 
system, represented by vectors y(i), might not synchro- 
nize because of instabilities on the synchronization man- 
ifold x(t) « y(i)- They introduced an approach that 
looked at many of the ideas of observers from the point 
of view of this instability. They also recognized that the 
Kalman-like filters that are often lifted from their natu- 
ral setting in linear problems to help achieve observability 
or synchronization in nonlinear problems are directed to- 
wards the reduction of noise in the tasks at hand rather 
than directed towards instabilities in the formulation of 
the problem itself. 

We have introduced a variational formulation for state 
and parameter estimation in nonlinear problems that fo- 
cuses on these issues as aspects of one problem. To il- 
lustrate the issues involved and the solutions we suggest, 
we have constructed an experiment on a simple, widely 
utilized, three variable nonlinear circuit, and we present 
both the formulation of our state and parameter estima- 
tion in that context as well as presenting the results of 
our analysis of the observed data for the circuit. 

The circuit is shown in Fig. [1] It is a damped LRC 
resonator together with a bipolar junction transistor pro- 
viding a nonlinear dynamical element. There are three 
independent dynamical variables for its description, and 
we choose them to be VcE{t), the voltage at the collector 
relative to the emitter, Vsit) the voltage at the emitter 
relative to ground, and the current through the in- 
ductor. 

Our challenge in this paper is to use measurements 
of VE{tn) taken every tn = nAt; At — 10/iS for n = 
0, 1, TV to estimate the parameters in the model of the 
circuit and to estimate the unobserved state variables 
ycE{tn) and/L(t„) over the same selection of observation 
times. Once we have estimated ^/^^(iAr), /lC^at) and the 
fixed parameters, we may use the model for prediction 
for t > tM- 

It is the transistor in particular that makes this an in- 
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FIG. 1: Colpitts circuit 



teresting system for study. A simple transistor model is 
well known, but what is perhaps not very well known is 
that real transistors will follow this model only approxi- 
mately. As we will show, an imperfect model can have a 
significant impact on our analysis, especially on efforts to 
extract accurate parameter values. We attempt to find 
methods to test the accuracy of our results and poten- 
tial improvements to the model. In carefully examining 
the effect of coupling strength between model and data 
we find interesting and anomalous behavior which may 
provide a useful test for the quality of different models 
for the system. 

We examine two methods for estimating parameters 
that build on the idea of stabilizing the synchronization 
manifold between the observations and the model dy- 
namics [l^l . 

• One method (initial value) determines the param- 
eters and the state variables at the initial time 
Ve(0), Fc£;(0), /l(0) by solving the model equa- 
tions for a selection of parameters and initial con- 
ditions to produce a model counterpart of the ob- 
served variable's time series. A metric for the dis- 
tance between the model output and the observa- 
tions is then minimized using a standard procedure 
in initial condition, parameter space. The advan- 
tages of this method include speed and flexibility 
of use. Software implementing this method is avail- 
able from one of the authors (PB) [25|. 

• The other method (constrained optimization) 
works in the larger space of model variables at each 
observation time i„, a coupling/control parameter 
at each of these times, and the time independent 
parameters. This method offers more efficient and 
precise control of synchronization, though it suffers 
from a number of practical problems, e.g. setup, 
speed, memory and numerical problems. Software 
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for implementing this method is also available [2£ 



A. Circuit Equations 

We work with a nonlinear oscillator of Colpitts vari- 
ety. This has a long history in the study and utilization 
of nonlinear oscillators for a variety of technological ap- 
plications. While it is likely that many engineers working 
with Colpitts and other nonlinear oscillators encountered 
chaotic oscillations, Kennedy [l^l appears to be the first 
to have recognized this was not noise or an unwanted 
defect in circuit behavior. We utilize a chaotic Colpitts 
oscillator for illustrating our suggestions on estimating 
parameters and unobserved states in a chaotic system. 
It is important to move beyond the numerical study of 
these oscillators and demonstrate in an experimental set- 
ting how the ideas work, because theoretical models at 
best only approximately describe real experimental sys- 
tems. 

We built the circuit shown in Fig. [T] using the the 
following components: Ci = 7.44 ^F, C2 = 7.23 ^F, 
L = 11.74 mH, Ree — 392 ft, and the power supply 
Vcc — 5.03 V, Vee — —5.10 V. The capacitors were non- 
electrolytic with a maximum voltage rating of at least 
30V. The inductor was air core to avoid any nonlinearity 
due to core saturation or hysteresis. We used a 100 
potentiometer in the circuit to adjust the parameter R. 
Note that R in all the equations includes the resistance 
of the potentiometer plus the resistance of the inductor, 
which was not negligible in our case. We used a 2N2222 
BJT small signal transistor. The fundamental frequency 
of this oscillator is approximately /o = 1 / {2Tr y/ LCeq) , 
where 1/Ceq = 1/Ci -I- I/C2. For the circuit elements we 
used, fo « 770 Hz. 

The dynamics is described by three coupled first order 
differential equations, obtained directly from the circuit 
using Kirchoff's laws: 



dVcEjt) 
dt 



hit) - IciVEit)), 



(1) 



dt 



Re 



lB(VE{t)), (2) 



L^^ = Vec~VEit)~VcEit)-RlLit), (3) 
dt 

where Ic{Ve) and Ib{Ve) are the currents into the col- 
lector and base of the transistor. The state of the circuit 
is completely described by the three dynamical variables: 
Vc E {t) the potential at the collector relative to the emit- 
ter, Ve(0 the potential at emitter relative to ground, and 
I Lit) the current through the inductor. A key ingredient 
needed to complete the description is the set of equa- 
tions that specify the transistor currents. Here we used 
a simplified version of the Ebers-MoU equations [H, [2^ , 



Ic{Ve) = /^exp 
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where (3f is the forward current gain and Is is the re- 
verse saturation current, both properties of the particu- 
lar transistor. Vt = kT/e is the thermal voltage. These 
equations are nonlinear and are what leads to the inter- 
esting behavior of the circuit. Inaccuracy in the transis- 
tor model has important consequences for our analysis. 
Later in this paper we discuss some ways to improve it. 

Although in this test circuit we can take experimental 
measurements of all three of the variables, it will often 
be the case in experiments on other systems that only 
some of the variables can be measured; this is certainly 
true, for example, of weather forecasting and nervous sys- 
tems. So in our analysis we will assume that only VE(t) 
is available for analysis, and that the circuit parameters 
are unknown. 



B. Circuit Dynamics 

The values of the circuit parameters determine the be- 
havior of the circuit. Although all the parameters mat- 
ter, here we focus on changing only the resistance R by 
adjusting the potentiometer, and keeping all the other 
parameters fixed. 

The simplest possible behavior of the circuit occurs 
when the three dynamical variables are constant in time; 
namely, a fixed point. To find this point, set the time 
derivatives in the equations of motion to zero. The emit- 
ter current, Ie = Ic + Ib, can be found from Eqs. Q 
and ([5]), and set equal to the current through the resistor 
Ree, 



Ie^Is{1 



exp(- 
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The fixed point emitter voltage V^'^ can then be found 
numerically for specific parameters. The other two fixed 
point values V^?^ and I^'' can then be calculated, using 
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Is exp(-Vi"Vl^T), and V^"j^ = Kc - RI, 
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The fixed point is a useful thing to calculate, because 
it can be used as a first step for comparing the model 
predictions to the voltages measured in the actual circuit. 
This fixed point is always a solution to the equations of 
motion, but it stable only when R is large enough. 

As R is decreased the fixed point becomes unstable and 
the circuit oscillates at approximately the fundamental 
frequency fa. As R is decreased further there is a series 
of period doubling bifurcations leading to chaos. This is 
shown in the bifurcation diagram Fig. O To construct 
this diagram, we integrated the circuit equations using a 
range of values for R, and plotted the value of Vce only 
at the times t„ where both VE(t„) = Vth and '^^"^f"'' > 0. 
The threshold value of Vth — —0.6 V was chosen because 
that is about where the transistor switches from "on" to 
"off'. 
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FIG. 2: Bifurcation diagram calculated by numerically inte- 
grating the model equations over a range of values for R. The 
value of VcE{tn) is plotted for all tn where VE{tn) — -0.6 V 
with Ve increasing. There is a period doubling bifurcation at 
R X 59S1 and again at i? « 48f2. 



We recorded all three dynamical variables 
VcEit),VE{t), and /^(t) when the circuit was op- 
erating in a chaotic regime. The orbit traces out a 
strange attractor in this three dimensional space. A 
two dimensional projection of the attractor is shown in 
Fig. [3] with data both from our experiments and from 
our model. 



II. PARAMETER AND STATE ESTIMATION 
USING SYNCHRONIZATION 

A. Background 

Although we were able to measure all the dynamical 
variables for our Colpitts oscillator, we use this as an 
example system to evaluate methods for estimating un- 
known parameters and unobserved state variables. We 
proceed as though we know many of the 'external' pa- 
rameters driving the circuit, such as the constant voltages 
Vee and Vcc while treating the parameters within the cir- 
cuit, especially those of the model used for the nonlinear 
operation of the transistor, as unknown. 

We proceed, then, with the assumption that we observe 
only VE{tn) where t„ = nAt, with the data time-step 
At = 10/iS. We wish to use this information to evalu- 
ate the remaining parameters and the value of the state 
variables Vcsitn) VE{tn) and Ihiin) for n = through 
N . By assumption we have data for VE{tn), but this will 
be slightly different than the output from the model due 
to noise and imperfections in the model. Except where 
otherwise specified, our results are for the case N = 1000 
or t^r = 10 ms. For all of the analyses in this paper, the 
model uses both a half-step and a full-step for improved 
accuracy. Often the model half-step h will be chosen to 
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FIG. 3; Two dimensional view of the attractor: Top Panel 
Vsit) versus lL{t) calculated by integrating the model equa- 
tions; Bottom Panel VE{t) versus I hit) as observed in our 
experiments. 



be equal to the data time-step. i.e. h = At = 10/iS. 
In this case the full-step will be twice this value and the 
number of full-steps used in the analysis will be N/2. 
As is discussed later in this paper it is sometimes ad- 
vantageous to choose a smaller value for /i, typically by 
selecting a value for which some multiple of h is equal to 
At. 

It is very important to choose an appropriate quantity 
of data to analyze. With too little data there may be in- 
sufficient information to determine the unknown parame- 
ters, while with too much data, not only is the processing 
slowed down, but (as we will show) the problem of main- 
taining synchronization becomes more severe. Once the 
fixed parameters are estimated along with VcEitN) and 
lL(tN), we may use them and the observed VE{tN) as 
initial conditions for the three dynamical circuit equa- 
tions to predict VE{t > tpf). We would also predict 
VcE{t > ti^) and /^(i > tN) at the same time, but in 
the scenario we are envisioning, those are not observable, 
so we would have no way to verify the predictions. The 
advantage of our experimental setup and our measure- 
ments of all state variables is that in the case we are 
discussing here, we do, in fact, know VcEitn) and /^(in) 
and can further examine the validity of our estimates. 
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In the interesting physical settings we have in mind, this 
will not be the case, and we wish to have some sense of 
the reliability of the methods in that instance. 

We make use of the phenomenon of synchronization 
of chaotic systems. The key idea is that two coupled 
systems will synchronize the most readily when they are 
identical. In our case, one system is the "driver," which 
here is the actual circuit, and the other system is the 
"response" or model system which is implemented nu- 
merically using the model equations with an additional 
term to couple to the driver system. The response system 
has adjustable parameters, and our goal is to find the val- 
ues for the parameters that best match the parameters 
of the driver system. Here we introduce a simple form 
of coupling by taking the difference between a measured 
data variable and the corresponding model variable at 
each time step, multiplying this by a coupling strength 
u, and adding this expression to the right hand side of the 
corresponding differential equation. This can be done for 
as many variables for which we have data, but since we 
are using Colpitts as a "test system" we choose to only 
apply it to Ve ■ This results in the following sets of equa- 
tions, in which the primed variables represent our model 
of the system (response) and the unprimed variable Ve 
represents the experimental data (driver): 

dt C[ 

^ ^ + 

^ -Ve- V6e - R'I'l 

Another possible coupling scheme (that we have not 
used) would be to replace V^{t) with aVE{t) -I- (1 — 
a)V^(<) on the right hand side of all three circuit equa- 
tions, where < a < 1. 

When model and data are different at some time-step, 
i.e. when VE{t) ^ V^{t), the coupHng term will try to 
move in the direction that closes this gap. Note also 
that the coupling strength u can be allowed to change 
with time-step or it can be independent of time-step. It 
can also be put under the control of the optimization 
program, or it can be fixed, or it can be manually changed 
through a set of values. We explore all of these options. 

Note that numerical instability is likely if u is allowed 
to be too large. The limiting value is approximately l/2h, 
i.e. one should require 

u < l/2h. (8) 

To increase u beyond this point one can make a corre- 
sponding reduction in the model half-step, possibly re- 
quiring the experimental data to be interpolated. An- 
other common reason for reducing the timestep would 
be when it is necessary to improve the accuracy of the 
numerical integration of the equations when no coupling 
is involved. 



An alternate scheme that avoids the need for interpola- 
tion of the data and that we have used quite successfully, 
is to remove the synchronization term from the differ- 
ential equations, and instead make small synchronizing 
corrections directly to the model variable. These correc- 
tions are applied between integration steps when there is 
data available at that time point, discontinuously moving 
the model variable some fraction of the way towards the 
corresponding experimental data value. The small steps 
in the model variable will tend to vanish on approach to 
perfect synchronization. This method of coupling can be 
shown to be essentially equivalent to the other method if 
that fraction is 1 — exp(— wAt). Using subscripts — and 
+ to represent before and after the discontinuous step, 
we have for our current problem: 

Fi;(t„+) = V^itn-) + exp{~uAt))iVE{tn) - ■(4(i„-))(9) 

Note that this removes the stability problem for 
large u, i.e. Eq. ([8]) no longer applies. The case u 
goes to infinity corresponds to a jump in the model vari- 
able all the way to the data value for that time-step. One 
way to obtain Eq. ([9]) is to leave the coupling term in the 
differential equation but replace u with uAt J2 ^(^ — ^n), 
where the sum is over all data points (indexed by n). At 
is the data time-step and (5() is the Dirac delta function. 
This can be thought of as a special kind of time depen- 
dent coupling. The delta function can be integrated out 
to give the fractional displacement rule above. It is of- 
ten the case that the desired timestep for the model is 
smaller than what is available or possible to achieve for 
the experimental data. In cases like this where the data 
is sparse, the method just described allows the time- 
step of the data to be a multiple of the time-step of the 
model, conveniently solving this problem without need- 
ing to interpolate the experimental data and avoiding the 
additional inaccuracy that would be introduced by that 
interpolation. 

We require a metric in which to compare our observed 
VE{tn) = VE{n) = xi{n) with the estimates we generate 
from solving the model equations for V^(tn) = V^{n) = 
yi{n). These estimates depend on the parameters, which 
we call p now, and the initial conditions y(0) = {Ve(0) = 
2/i(0),V'ci;(0) = 2/2(0), /i(0) - y3(0)}. We adopt the 
usual least squares or minimum variance metric as the 
distance between the experimental observations and the 
value of our model output 

1 ^ 

C(p,y(0)) = — -y^(a:i(t„)-yi(t„,p,y(0)))^ (10) 

We call this the cost function or objective function, and 
it is a measure of the quality of synchronization. 

The goal is to find the parameters and initial condi- 
tions which minimize the cost function. One problem is 
that the path through parameter space from the initial 
guess to the best result may contain local minima which 
will prematurely end the minimization process resulting 
in an incorrect result. But even when this is not the case, 
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FIG. 4: Cost (Eq. (fTTjl ) as a function of R' for four different 
coupling strengths (decreasing from bottom to top). The cal- 
culation was done over 10 ms, with the two systems starting 
at the same point in phase space, and all parameters identical 
except R^ R! . Increasing u smoothes out the surface, but it 
also makes the minimum less well defined as u becomes large. 



the "downhill" path through a high dimensional param- 
eter space can sometimes be extremely convoluted and 
lengthy. We can see an indication of this in Fig. [1] by 
looking at the variation of the cost function, C(i?'), 

1 ^ 

when only the resistance R' in the response circuit Eq. ([7|) 
is varied with all other parameters held fixed. The curve 
labeled u = corresponds to integrating the model equa- 
tions and evaluating the cost function Eq. (fTO|) . The 
other elements of Fig. [4] will be addressed in a moment. 

This complicated behavior of C{R') impedes one's 
ability to seek the global minimum where we expect 
that R' will equal R of the driver circuit-the source 
of our data. As one might expect, as more param- 
eters are added, the graph becomes a very complex 
hyper-surface that we must explore. This complexity 
is associated with the instability of the synchroniza- 
tion manifold x(i) = [^bW, ^CE W, W] = = 
^c_e(*)j (0] ill the six dimensional state space 

[x(0,yW] = [VE{i),VcE{i)jL{t),Vi,{t),Vl;E{t)J'L{t)] 

of the coupled oscillators. This instability leads to ex- 
treme sensitivity of the orbit y(i) on its parameters, in- 
cluding here i?', and, as it happens on its initial condi- 
tions y(0) though the latter is not shown here [lOl- This 
complexity of the surface of a cost function in parameter 
space has been known for some years [31, 32]. 

It should be noted that when additional parameters 
are included, many (perhaps most) of the local minima 
shown in Fig. 2] will no longer be minima. That is to 
say, in many cases it may be possible to escape from 
one of these apparent minima by simply moving in some 



orthogonal direction. It also presents a somewhat over- 
simplified picture, in that the real experimental data will 
often include some significant noise, and also the model 
we are using to represent the experimental system may 
have significant imperfections. As a result no choice of 
parameters will produce perfect synchronization and the 
optimal cost will be nonzero. 

With these considerations in mind, now we examine 
the other curves in Fig. 3) They correspond to increas- 
ing values of u, and we see two effects as u is made larger: 
(1) the complex, rough nature of C{R') as a function of 
R' becomes smooth, and (2) the magnitude of C(i?') de- 
creases. The origin of these two effects comes from the 
synchronization of the data xi{n) — s- yi{n,R'). As u 
becomes very large, it drives xi(n) — > yi{n,R') approxi- 
mately as i leading to C{R') decreasing as Interest- 
ingly, for small or moderate size u, there is structure in 
the dependence of the cost function on u that can reveal 
information about the quality of the model, the presence 
of noise, and other attributes useful in interpreting the 
outcome of the numerical optimization. 

The smoothing out of the surfaces in parameter and 
initial condition space is attributable to the reduc- 
tion of any positive conditional Lyapunov exponents 
(CLEs) [11, Si,! by the added term u{VE{t) - V^{t)). 
Note that the term "conditional" is used because the ex- 
ponents depend on the experimental variable Ve which 
is treated as an external drive signal whose dynamics 
are not explicity considered. The CLEs are evaluated by 
iterating the Jacobian DF(y(<)) along an orbits of the 
dynamical equations, and this change in the dynamics 
takes the original Jacobian and modifies it by subtract- 
ing u from the (2,2) component in our example. This 
allows us to reduce the one positive CLE to a negative 
value, synchronize the data and the model output, and 
make the contributions to the cost function coherent. 

From the point of view of the nonlinear dynamics ex- 
pressed in Eq. ([7]), we are, by coupling the equations to 
the data, reducing the largest conditional Lyapunov ex- 
ponent from positive to negative. This idea is illustrated 
in Fig. O which shows how the largest CLE depends on 
the coupling strength, and depends on which variable is 
being coupled. Coupling to /^(t) will not cause synchro- 
nization, at least for m < 1. 

It is desirable to have m — > at the end of any es- 
timation procedure as the term u{VE{t) — V^(t)) in our 
equations is not a result of Kirchoff 's laws or any physics, 
but is introduced to regulate instabilities on the synchro- 
nization manifold. To use the equations for prediction, 
it should be gone from the dynamics. 

Although increasing u appears to make our parameter 
finding problem much easier it can also have an undesir- 
able effect. This comes from the possibility that in real 
experimental systems, where the model is imperfect, the 
coupling term in the equations of motion will not com- 
pletely vanish for any choice of model parameters. This 
in turn can cause the parameter values that minimize the 
cost function to be dependent on the magnitude of the 
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FIG. 5: Largest conditional Lyapunov exponent (CLE) as a 
function of coupling strength. With coupling to Ve or Vce 
the largest CLE is reduced from positive to negative, but for 
coupling to II the largest CLE is actually increased. 



coupling. To deal with this, we can initially use a very 
strong coupling to arrive at a first approximation to the 
desired parameters, and then the coupling can be reduced 
gradually to smaller and smaller values, either to zero or 
to a value where its impact is small. To accomplish this 
it is sometimes efficient to allow the optimization process 
to control the value of u. In this case it may sometimes be 
necessary (though not always) to add a "coupling cost" 
to the cost function Eq. (fTT|) so that the cost can decrease 
as the coupling is decreased. 



C{p,u) = 



N- 



N 
n=0 



ixi{tn) - yi(in,p) 



2 2 
77 U , 



(12) 

where 77 is a number assigning a weight to the synchro- 
nization term relative to the coupling term, and p is a 
vector comprising all of the parameters and initial condi- 
tions. In some cases we allow the coupling to depend on 
time-step, based on the idea that the strength of the cou- 
pling required to make synchronization occur may vary 
across the system attractor. One might expect that in 
certain locations the data variable being measured may 
sometimes be in good alignment with direction of max- 
imal stretching of the attractor. This alignment can be 
examined explicitly through the Lyapunov direction vec- 
tors (in Rcf. [23] see discussion starting on page 2796). 
This does of course complicate the optimization problem; 
if there are 500 full time-steps for the model, and we wish 
to optimize this u(t) then we have just added an addi- 
tional 500 unknown parameters to the problem. We give 
results for both cases. 



B. Effect of Positive Lyapunov Exponents 

The Lyapunov exponents, particularly the most posi- 
tive one, can be very useful in determining how to analyze 
the experimental data. These measure the average rates 
of expansion or contraction of a chaotic attractor along 
different directions of the phase space. There are meth- 
ods for calculating these from experimental data [s^ [13] , 
and other, generally easier, methods that can be used 
when a good model is available and the parameters are 
known (in Ref. [IHl see discussion of QR decomposition 
method beginning on page 650 - software implement- 
ing this method is available [11]). When only the most 
positive exponent is needed, the problem is much easier 
- nearby initial conditions (on the attractor) will sep- 
arate on average as e'^* , where A is the largest Lya- 
punov exponent. For our system and parameter choices, 
we determined the most positive exponent to be about 
0.35(l/ms). For our choice of time interval tN = 10 ms 
for most of our results, one can expect that a small error 
in the initial condition at i = will grow by a factor 
e^-^ = 33 by the end of the interval. As a result of this 
moderate growth over our interval, it is relatively easy to 
reduce the coupling to zero without degradation to the 
fit to the experimental data. Going much beyond 10 ms, 
however will change this picture. In order to reduce the 
coupling all the way to zero we need to require that 



Mn < r. 



(13) 



where F is probably not a precise constant, but is rather a 
number that is roughly about 10 and which may depend 
somewhat on the details of the system and how hard 
we are willing to work on the analysis. An increase to 
100 ms definitely puts us over this limit for the present 
problem, drastically increasing the growth factor to e'^^ = 
1.6 X 10^^. In this case it is quite impossible to reduce 
the coupling all the way to zero without major loss of 
synchronization. We must choose to either reduce the 
time interval, or accept the problems that may possibly 
be caused by maintaining a nonzero coupling parameter. 
In the latter case, having a time dependent u{t) presents 
a clear benefit in being able to apply the coupling more 
efficiently at points where it is most beneficial. 



C. Optimization Methods 

We used two different optimization methods, an ini- 
tial value method and a constrained method, which we 
describe below. In both cases the primary goal was to 
find the best values for the parameters and all of the 
state variables based on the available experimental data 
which in this case is a time series for VE{t). Once this 
has been achieved we can also make predictions by inte- 
grating the model into the future, beyond the last data 
point. It turns out that for both of these methods, nu- 
merical accuracy of the model is improved by using a 
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half time-step, at which equations are evaluated, while 
the full accuracy is achieved only at the full time-step 
points. It was convenient for our data to use a half-step 
for the model that was equal to the time-step for the 
data. So we use i„ to represent the time points for the 
data as well as these half-step time points for the model. 
The even values of n correspond to the full-step time 
points. However, for other experiments it will often be 
the case that this will not be a good choice. In general, 
the model half-step should be chosen to be as large as 
possible without adversely effecting numerical accuracy. 
But, often there are limitations on the sampling rate for 
experimental data that require that the data time-step 
be larger than this value, although usually one can ar- 
range for the data time-step to be some multiple of the 
model half-step. 

The first method we examined is an "initial value 
method" similar to the one discussed in Voss et al. [33|, 
but uses a time independent coupling term to synchro- 
nize the experimental data set to the model. This method 
treats the initial conditions y(0) as additional parameters 
to be determined. The cost function can be computed at 
any point in this extended parameter space by numeri- 
cally integrating the equations of motion with those val- 
ues of the parameters, and comparing the result to the 
data using the cost function Eq. (fTU)) . We use standard 
minimization methods to search the parameter space for 
the minimum of the cost function. Integration of the 
equations is typically done using the standard fourth or- 
der Runge-Kutta algorithm. Note that this algorithm, 
as mentioned previously, involves function evaluations at 
the midpoint of each full time-step. Software (DataSync) 
implementing this initial value method is available from 
one of the authors (PB) on request [25j . 

Two of the minimization methods that we tried with 
the initial value method were the Broyden-Fletcher- 
Goldfarb-Shanno (BFGS) "quasi-Newton method" [13, 
and the Brent-Zangwill-Powell-Smith (BZPS) "direction 
set method" [s^. We first set upper and lower bounds 
and starting values for each of the parameters to be used 
by the minimization method. We used physically rea- 
sonable values for bounds and starting values. Most of 
the results in this paper were obtained using the BZPS 
method. BZPS is convenient because it does not require 
evaluation of derivatives of the cost function. The deriva- 
tives required by BFGS can however be determined nu- 
merically. It has been our experience that BZPS is the 
most reliable method, even though it may sometimes 
be shghtly slower than the BFGS method. The BZPS 
method works by performing line searches to find the 
minima of the cost function along specific directions in 
parameter space. Each direction is searched along in a 
specific sequence, always staring off at the location in pa- 
rameter space where the last search ended. Initially the 
search directions are just the cardinal directions, but as 
the process goes on the directions are updated based on 
the results of previous iterations. 

In this initial value method, we often set the coupling 



strength u to be time independent and fixed at a user 
specified value. We use the method in several steps. Ini- 
tially we use a large u to get from the initial guess of 
parameters to the correct region of the parameter space. 
Then we repeat with smaller u, starting where we left off 
in parameter space, in order to refine the search. This 
approach is suggested by Fig. 01 which shows the that 
the cost function becomes smoother with large coupling. 
Provided the time interval is not too long, as discussed 
previously, u can be reduced all the way to zero without 
significantly impacting the fit between model and data. 
When this is not the case, there will be an (approximate) 
lower limit to u below which synchronization is lost. It is 
also possible to include u in the cost function (Eq. (fT^ 'l. 
and have it be reduced as part of the optimization pro- 
cess, although, as will be shown, this must be set up 
carefully in order to get the best result. 

The second approach we have utilized for this nu- 
merical optimization problem we call the 'constrained 
method'. In this method, one treats the values of the 
model state variables y(tn) for all even n on an equal foot- 
ing with the parameters p. The differential equations are 
discretized to get equality constraints relating y(tn) to 
y(^ji+2), effectively forcing the optimization algorithm to 
find values which satisfy those differential equations. The 
discretization process involves evaluation of the equations 
at the midpoint of the full time-step, i.e. the points t„ 
where n is odd (the details of this "Hermite-Simpson" 
process can be found in Ref. [IBl)- As before, we were 
able, for our particular data set, to make the choice that 
the half time-step for the model was equal to time-step 
for the data. As a result, when we are using A^-|- 1 exper- 
imental data points in our analysis, this method needs 
to minimize the cost function in a space of dimension 
D{1 -\- N/2) -t- L, where D is the number of variables and 
L is the number of parameters. When we allow the cou- 
pling u to depend on this will add another 1 -I- N/2 to 
this total. Thus this dimension is typically about 2000 
when we use iV + 1 = 1001 data points, i.e. the space 
we are optimizing over is extremely large. We have suc- 
cessfully used up to -I- 1 = 10001 data points in the cal- 
culation. We implemented the constrained method using 
the constrained nonlinear optimization software package 
SNOPT with the optimal control interface SNCTRL 3^1 . 
This is available publicly [l^. We have used a symbolic 
mathematics code (Matlab or Mathematica) to provide 
the derivatives of all equality constraints and of the cost 
function. 

The output of the numerical optimization using ei- 
ther method is an estimate of all the parameters along 
with an estimate of all of the unobserved state variables 
2/2 (in), 2/3 (in ),■•■, J/D (in) ovcr the time interval of data 
presented to the model. For the Colpitis experiments, 
this gave us an estimate of the parameters not held fixed 
(see below) as well as VcE{tn), Ibitn) ; n = 0, 1, 2, N 
over the interval from to = to — NAt. From the 
values of VE(tAr), VcE{tN), IhiS-N) and the parameters we 
can use the differential equations to predict the temporal 
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development of the oscillator. 

As stated previously, we often wish to reduce the cou- 
pling to zero, at which point the rms error in fitting the 
data will be a good measure of our success in finding the 
correct parameters. However, if it is impossible to reduce 
the coupling to zero, perhaps because of the excessive 
length of the data set, we can use the information from 
either method to ask an important consistency question 
about our estimation results by examining the magnitude 
of the ratio 



FE{t? 



FE(tY + [u{t){VE{t)-V^{tW 



(14) 



where Fe(<) = dV^/dt without the coupling term (see 
Eq. ([7])). If R{t) is near unity over the time segment 
where we estimate parameters and state variables, the 
estimation is consistent as the coupling term is small 
compared the the dynamics in FE{t). If, however, R{t) 
departs from unity, then the accuracy of any synchro- 
nization x(i) Ri y{t) is due to the coupling. If R{t) w 1, 
we call the model consistent with the data. If not, espe- 
cially if R{t) is near zero, we call the model inconsistent 
with the data. 




2 4 6 

Coupling Strength [1/ms] 

FIG. 6: Optimized rms error as a function of fixed, time in- 
dependent coupling strengtli using 10 ms of real data and 
analyzed with the standard model. The result is anomalous - 
the initial increase in error with increasing coupling was un- 
expected. Note that the coupling strength is not dependent 
on time-step for these results but rather is fixed at the values 
given on the horizontal axis. Note also that the parameters 
are re-optimized for each point on the curves to reduce the 
rms error as much as is possible. 



III. STATE AND PARAMETER ESTIMATION 
USING THE CHAOTIC COLPITTS CIRCUIT AS 
A TEST SYSTEM 

As we have described, we operated the Colpitis cir- 
cuit in a setting where chaotic trajectories for the states 
as a function of time were generated. We recorded 
VcE{tn),VE{tn), Ihitn) at intervals of At = 10 ^s for 
100 ms {tn = nAt,n = 0, 1, 10"*). We then (un- 
less stated otherwise) used the first 10 ms of the sam- 
pled VE{tn) data as input into our estimation proce- 
dures. As indicated earlier, we held fixed several pa- 
rameters: Ci^VccVee, and i?ee- There remained six pa- 
rameters to determine from the 10 ms of VE{tn) data: 
C2,L, R, Is, Vt, Pf- Instead of directly treating the very 
small quantity 1^ as a parameter, we defined Vq = 
—Vt log tJ, where /o = 1 mA, and searched on Vq- 

We measured all values of the circuit elements directly, 
and we found the transistor parameters by measuring the 
I-V curves, /c(^B£;) and Ib{Vbe), so that we could com- 
pare to the parameter estimates we obtain. The resolu- 
tion of the VE{t) measurement was AV — 0.01 V. If 
the only source of uncertainty in the experiment was due 
to resolution of the measurement, we would expect an 
rms error (square root of the cost function, Eq. ((T0| ) of 
yJC = AV/^/vi « 0.003 V. In both methods described 
below, the actual rms error between the model and the 
data comes out to be about four times larger, suggesting 
there are other sources of error. 



A. With Time Independent Coupling 

Here we use the DataSync software to implement the 
initial value method with time independent coupling u. 
The response system is now a model of the circuit with 
six unknown parameters plus the three unknown initial 
conditions Ve(0), VcbCO), and /l(0). Exploring the be- 
havior of the rms error in fitting the data as a function of 
coupling strength u and find a very anomalous result 
shown in Fig. [51 We expected the error to decrease with 
increasing coupling strength. This is in fact the case for 
u > 1.7. However, for u between and 1.7 the opposite is 
true - the error rises and with a very significant positive 
slope even near u — 0. 

Recall that when we wish to have the optimization pro- 
cess reduce the value of the coupling u towards zero, we 
would normally include a term containing the coupling 
strength in the cost function as in Eq. However, 
if the coupling is restricted to the anomalous range be- 
tween zero and the value where the rms error starts to 
decrease, then the coupling term may be omitted from 
the cost function entirely. We tested this using the initial 
value method with the coupling restricted to the range 
from 0.0 to 1.0 (1/ms), and with all parameters set ini- 
tially to values that were either 1.4 or 0.7 times their 
best values and with the coupling started at 1.0. The 
optimization was successful and the coupling ended up 
at about 3 x 10"^. 

If a wider range of coupling is desired then the value of 
rj that appears in the cost function must be chosen large 
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enough to counteract the downward slope that occurs for 
large coupling. If we let e{u) represent the optimal rms 
error as a function of u, then the optimal cost function is 
given by C(u) = e('u)^ + rj^u'^. Note that optimal means 
with parameters adjusted to minimize the cost for the 
particular choice of u. If this function has a minimum for 
positive u then the optimization process could get stuck 
there. To look for a minimum we set the derivative equal 
to zero resulting in the equation: 7]^ — {~de/du){e/u). 
We can eliminate all such minima if we choose 77^ so that 

ri^ > i~de/du){e/u), (15) 

for all u in the desired range. In that case u should tend 
to go to zero. 

In the present case we can estimate de/du from Fig.[6l 
We found that for u in the range from 0.0 to 10.0 the 
minimum value for 77^ should be about 2.6 x 10^^. Again 
we tested this with the initial value method. With rj^ set 
to 2 X 10~^ the optimization got stuck in a minimum near 
M = 6, but when we increased rj^ to 3 x 10^® it no longer 
got stuck and successfully completed the optimization. 
For this run u ended up at about 1 x 10~^. 

Normally when starting to work on a new problem 
de/du will not be known in advance, so we would have 
to guess based on our knowledge of the problem and pre- 
vious analyses of similar problems. The difficulty in se- 
lecting T] can, of course, be avoided entirely by requiring 
u to take on a series of values ending with zero. 

As mentioned previously, if we wish to be able to re- 
duce the coupling to zero, we must ensure that the time 
interval t^r over which we are analyzing the data is not 
too large. In particular we must require that ijvA is not 
very large, where A is the largest positive Lyapunov ex- 
ponent (0.35 for our data). When t^q is 10 ms as it is 
for most of our results there appears to be no problem 
in reducing the coupling to zero. However if we increase 
to 60 ms this is no longer the case. When we try 
to reduce the coupling in that case we reach a value be- 
low which it becomes impossible to maintain a good fit 
to the data. This is illustrated in Fig. [T] To create the 
curve in this figure, we start by optimizing at a fairly 
large value of coupling, e.g. u = 10. We then reduce 
u in steps, each time starting the parameters with the 
values obtained in the previous step. At each step we 
re-optimize the parameters and make sure that the rms 
error is only slightly different from the previous step. If 
the step taken is too large the error settles on a value that 
is much larger than the previous step, often by a factor of 
10 or more. When we are near or below u = 1.0 this pro- 
cess starts becoming exceedingly difficult. We must take 
smaller and smaller steps and the optimization process 
becomes slower and slower. Soon it becomes impractical 
to make any more progress. For the case shown in Fig. [7] 
this occurs around u ~ 0.57, which is the lowest point 
shown. Thus for this data-set we cannot reduce 
the coupling to zero. Our choices are: (a) use the re- 
sults obtained at the minimal coupling u = 0.57, (b) redo 
the analysis using a smaller time interval so that we can 
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FIG. 7: Optimized rms error as a function of fixed, time in- 
dependent coupling strength using 60 ms of real data and 
analyzed with the standard model. Due to the length of the 
time series it is not possible to reduce the coupling strength 
to zero without a loss of synchronization. The points were 
found by moving downwards slowly in u and re-optimizing 
after each step. It is not possible to continue the curve below 
the marked end point - instead the optimization is forced to 
jump to a very unfavorable state with at least 10 times the 
rms error. 



reduce the coupling to zero, (c) redo the analysis using a 
time dependent u{t) as discussed in the next section, (d) 
break up the time interval into smaller subintervals, each 
of which is short enough to allow zero coupling. The last 
case would be classified as a "multiple shooting" method. 
Each subinterval would have its own initial conditions to 
be determined, so there would be an increase in the total 
number of unknown parameters. But this would be offset 
by the benefit of having more data points available 

The parameters calculated with time independent cou- 
pling and using the shorter time series of 10 ms are com- 
pared to the measured values in Table [H The parame- 
ter estimates mostly come out near the measured values, 
with the exception of the transistor parameters (3f and 
Vq. We explore this discrepancy later with an improved 
transistor model. 

We also compute a value for the uncertainty of the pa- 
rameters. This was calculated by doing exactly the same 
fitting procedure, with the same initial guesses for pa- 
rameter values, over nine different 10 ms segments of the 
same data set. Since all the segments of data came from 
the same circuit, we expect all the parameter estimates to 
come out almost the same, but the initial condition esti- 
mates to be different. The column labeled 'Uncertainty' 
in TableUis the standard deviation of the mean of each of 
the parameters over the nine different 10 ms segments of 
data. This is just a measure of the uncertainty from the 
numerical fitting process. There is an additional uncer- 
tainty presumably due to experimental noise and other 
effects which are not included in the model. 
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Name 


DataSync 


Uncertainty 


SNOPT 


Measured 


Units 


C2 


6.98 


0.03 


7.02 


7.23 




L 


12.28 


0.02 


12.2 


11.74 


[mH] 


R 


40.38 


0.01 


40.0 


39.3 


[Q] 


Vo 


0.663 


0.004 


0.661 


0.63 


[V] 


Vt 


25 


1 


25.0 


27 


[mV] 




74 


4 


72.0 


180 


[1] 



TABLE I; Parameters found by the two methods and the 
measured values. 



Name 


Measured 


Prediction 1 


Prediction 2 


Units 




-0.696 


-0.723 


-0.702 


[V] 




2.24 


2.31 


2.25 


[V] 


AO) 


11.12 


11.02 


11.16 


[mA] 



TABLE IL Fixed point values. The first column are the values 
measured directly from the circuit. The second column is the 
calculation of the fixed point using the fitted parameters using 
the standard model. The third column is the calculation of 
the fixed point using the fitted parameters from the improved 
model. 
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FIG. 8: Time dependent coupling u{t) for r; = .10 [V ms] 
(middle) and 77 = .14 [V ms] (bottom). The structure is very 
similar in both cases, but the magnitude scales approximately 
as u{t) (X r)~^ . When u{t) is allowed to be positive or negative, 
additional structure becomes visible(top). 



Once the parameters are estimated with the fitting pro- 
cedure, the model can be used to make predictions. One 
simple prediction is the location of the fixed point. The 
fixed point location is calculated using the estimated pa- 
rameters (but substituting the measured value of R, since 
the resistance was increased to make the fixed point sta- 
ble), and also can be directly measured from the circuit 
by making R large enough to make the fixed point stable. 
The measured and calculated values are shown in Table 
im The values predicted using the fitted parameters are 
within 5% of the measured values. The third column in 
the table contains predictions from an improved model 
that will be explain later in the paper. 



B. With Time Dependent Coupling 

Now we set the coupling u to a function of time 
u — > u{tn), SO the strength of the coupling (or control) 
to stabilize the synchronization manifold varies over the 
attractor. We use the SNOPT software as a way to 
implement the constrained method for solving the nu- 
merical optimization problem. As mentioned previously, 
for 10 ms of data we have space of approximately 2000 
dimensions to search. As SNOPT is used regularly for 
optimization in spaces of high dimension, this was well 
within its tested regime of validity. We presented the 
software the data Vsitn) at 1000 points, provided it with 
initial guesses and ranges for the six parameters, and set 
the initial guesses for all the state variables and coupling 
variables to zero. 

In the previous section we noted that when the cou- 
pling strength was below the peak in Fig. [S] there would 
be a tendency for it to move towards zero even with i] 



set to zero. For time dependent coupling we find that 
there are subintervals in which the coupling will tend to 
decrease to the lower boundary (usually zero) and others 
in which it will tend to increase often going to the upper 
boundary (recall that we must require that uAt < 1 for 
numerical stability). With a nonzero rj there emerges a 
set of peaks, shown in Fig. [8] whose amplitudes decrease 
with increasing rj. These peaks are not random, but are 
robust features for a given data-set. The optimal rms 
error e and the optimal cost C are now (in a contin- 
uous time picture) functionals of the coupling strength 
u{t). Optimal here means the value obtained when the 
cost has been minimized over the parameter space for the 
particular choice of u{t). For a time interval from to r, 
the cost can be expressed as: 

C[u] ^ e[u]^ + r uitfdt. (16) 
T Jo 

We can solve for u{t) by taking the functional derivative 
and setting it equal to zero: 

dC[u] 6e[u] W „ 

We define eo to be the value of e[u\ for u{t) — and 
we define Ko(t) to be the functional derivative of e[u] for 
u{t) = 0. For sufficiently small u{t) we may use these to 
solve for u{t): 

«(<)«-^^oW- (18) 
T 

Thus we find that the structure observed in u{t) is 
proportional to the functional derivative of the 
rms error with respect to the coupling strength. 
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The equation also shows that the amplitude of the 
observed structure in u{t) varies inversely with 77^ 

as can be seen in Fig.[8l We usually set the lower limit on 
u(t) to be zero, but when this is instead set to a negative 
value we can see additional structure corresponding to 
times when Ko{t) is positive. This is also shown in Fig. [51 

We found that u(t) is strongly influenced by noise and 
this is responsible for the jagged appearance of u{t). To 
explore this we generated u{t) using simulated data with 
added noise. We then repeated this with the same data 
but using a different noise sequence of the same rms am- 
plitude. This produced a dramatic change in the struc- 
ture demonstrating that the noise played a critical role. 
We then repeated the test again, this time using the orig- 
inal noise sequence but with double the rms amplitude. 
In that case we found that the structure remains essen- 
tially unchanged but the amplitude of that structure is 
increased by a factor of four, i.e. it is proportional to the 
square of the rms noise amplitude. 

It is also interesting to see what happens if we change 
the form of the coupling term in the cost function. In 
particular one can use a term that is linear such as ^u{t), 
where ^ is a constant, resulting in a cost functional of the 
form: 

CM =e[u]2 + i r \u{t)\dt. (19) 
T Jo 

If we take the functional derivative again we get 

Since this doesn't explicitly contain u(t) we can't set it 
to zero and solve as before. Note, however, that (in the 
small u{t) limit) the right hand side will be positive when 
^ > —2TeoKo{t). In these cases u{t) will go to zero (the 
presumed lower boundary for u{t)). For those values of 
t for which the opposite is true, u{t) will move upwards, 
either to the upper boundary, or to a value where our ap- 
proximation breaks down because u{t) is no longer small. 
Further, one may expect that there will be a critical value 
of ^ above which u{t) will be zero for all t. 

Note that in all of the cases above, we have been dis- 
cussing the optimized value of u{t). In the process of 
getting to that optimal solution, u{t) may temporarily 
take on much larger values as determined by the opti- 
mization algorithm. 

In addition to noise, the structure in Fig. [5] may be 
partially related to the varying relationship between the 
direction associated with the (scalar) data variable Ve 
and the Lyapunov direction vectors (in Ref. js^ see dis- 
cussion starting on page 2796), particularly direction of 
maximum stretching of the attractor that is associated 
with the largest Lyapunov exponent. At certain loca- 
tions on the attractor these directions are nearly parallel 
and so the coupling there is relatively effective in achiev- 
ing synchronization. This becomes increasingly impor- 
tant for long data-sets, particularly those which violate 



the condition of Eq. p^ . These cases are also ones 
for which there may be a significant advantage for 
time dependent coupling. We previously examined a 
case for which the time interval was 60 ms, compared 
to 10 ms for most of our results. In the results, shown 
in Fig. [71 we were unable to reduce the time indepen- 
dent coupling below 0.57. At this point, the rms error 
was 0.0292, considerably higher that the value 0.0146 ob- 
tained in a shorter data-set. Running this data-set with 
a time dependent u{t) we were able to reduce this to 
0.0161. The average u{t) is now much less than 0.57, but 
it contains a number of peaks with height of order 1, for 
which, at the time when they occur, the primary Lya- 
punov vector is relatively parallel to the coordinate axis 
associated with the measured experimental data variable. 
When using the linear form of the coupling term in the 
cost function, a similar result is found except that the 
few remaining peaks are very tall and extremely narrow. 
For ^ = 0.001 there are only 13 peaks in u{t) with an 
average height of about 6 and an average width of only 
slightly more than one full time-step. Everywhere else it 
is zero. One might reasonably assume from this, that for 
many problems with a good choice of ^ there will be little 
impact on the dynamics other than to maintain synchro- 
nization. Of course, if ^ is too large the synchronization 
will be adversely effected while if it is too small the num- 
ber of peaks in u(t) will become excessive and detrimen- 
tal as it attempts to correct for every noise bump in the 
data. 

In Fig. [51 we show that the peaks of the residual u{t) are 
localized in phase space. (Here we are using 1800 points 
in the data set and a quadratic term in the cost function) . 
We divided the phase space into cells and shade in each 
according to the average value of u{t) within each cell. 
We do the same thing a second time, but using VcE{t) as 
the data source instead of VE{t). The areas where u{t) 
is largest are now in different spots in phase space. 

In the cost function Eq. (|12p we explored various values 
of the weight 77 and report in Table U the values of the 
parameters for 77 = 0.1 [F][7ns]. If r] is chosen too small, 
u will become large and the coupling term will end up 
dominating the dynamics in Eq. This means that 
V^{t) will be forced to follow VE{t) even if the parameter 
values or the model itself is wrong. On the other hand, 
if rj is chosen too large, u will be driven toward zero and 
synchronization may not occur. In either extreme limit of 
T] very small or rj very large the ratio R{t) w 1, because 
either (Ve(^) — V^it)) or u will become very small in 
Eq. (HH). Therefore care must be taken to chose a value 
of rj that allows the two effects to be balanced. 

The parameter search methods provided consistency 
for estimates of the six parameters as displayed in TableHl 
In Fig.[TUlwe display both the observed values of the state 
variables [iLit), VE{t), VcE{t)] and their estimated values 
from the use of the constrained method. Only Ve(<) as 
observed was presented to the estimation algorithm, yet 
the estimation of the two unobserved state variables II (t) 
and VcE{t) is also quite accurate. 
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FIG. 9: Localization of peaks of u(t) in phase space. Phase 
space was divided into 40 x 40 x 1 ceUs ( 40 in the II and Ve 
directions, 1 in the Vce direction), and the average value of 
u{t) in each cell is indicated by the shading. The calculated 
phase space orbit is overlaid for reference. The coupling was 
either to Ve ( Top Panel), or to Vce (Bottom Panel). In 
both cases rj = .01 [V ms] and 1800 data points where used. 
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FIG. 10: Using 1000 measurements of VE{t) over 10 ms (dots 
in top panel), we estimate the other two state variables, 
VcEit) and /i(t) (solid lines). Also shown for comparison 
is the measurement of VcE(t) and lL{t) (dots in middle and 
bottom panels). 




FIG. 11; Using the first 10 ms of data for VE{t), sampled every 
10 fis, we estimated the parameters and the other state vari- 
ables VcE{t) and lL{t) over the interval < t < 10 ms. Then 
using the equations of motion with estimated parameters and 
the values of VE{t — 10ms) (as measured) and VcE{t = 10ms) 
and lL{t = 10 ms) (as estimated), we predicted the three state 
variables for 10 ms < t < 25 ms. The predictions (solid lines) 
are compared to the measured data (dots). Initially the pre- 
diction matches the measurement, but the two diverge after 
about 8 ms because the dynamics are chaotic. 



From the output of the estimation program we have 
an estimate of the values of the three state variables at 
t — 1000 At = 10ms, and using these as an initial con- 
dition for the model of the Colpitis circuit along with 
the values of the estimated parameters (Table IJ), we pre- 
dicted the behavior of the circuit for t > iiooo- These 
predictions are displayed in Fig.[TT] Each shows that the 
predictions from the estimated state of the oscillator are 
accurate for about 8 ms or so beyond the point in time 
t = 10 ms where the state estimate was made. This is the 
kind of prediction one finds when a dynamical system is 
chaotic as small errors in parameters or initial conditions 
are exponentially magnified in the behavior of the orbit 
because of the positive Lyapunov exponents of the sys- 
tem. What is important here is that we have been able, 
through our estimation methods, to accurately enough 
evaluate the state variables that were not observed that 
our model predictions of the experimental system show 
high accuracy. 

Because we use a coupling to the synchronization man- 
ifold as part of our estimation procedure we need to ask 
whether the synchronization we observe is the result of 
large coupling. At the end of the optimization process 
we calculate the ratio R^{t) from Eq. fH]) . and see that 
R^(t) > .99 at all times as shown in Fig. [121 This means 
that the coupling has a negligible impact on the dynam- 
ics, as it should, and that our model is consistent with 
the data. 
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TABLE III: Parameters found with the improved transistor 
model, which includes emitter resistance Re, compared to the 
measured values 
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FIG. 12: The consistency check R^{t) defined in Eq. (|T4)l . 
Since R? (t) « 1 the model is consistent with the data. In this 
calculation rj — .1 [V ms]. 



C. Effects of Model Deficiencies and Results for 
tlie Gummel-Poon Model 

The Ebers-MoU model is only an approximate model 
of a real transistor. This causes interesting challenges for 
the parameter calculation process and makes the Col- 
pitts oscillator more interesting as a test problem. In a 
general problem, it may often be the case that the value 
of the parameters themselves will be of interest. But 
in the process of trying to optimize the output of an 
inexact model sometimes the deficiencies in the model 
may be partially compensated for with parameter val- 
ues that are dramatically shifted away from their true 
values. In the current problem this may have occurred 
with the current gain parameter (3f, for which we ob- 
tained a value of about 72 by minimizing the cost func- 
tion compared to a measured value of about 180. An- 
other anomaly was observed in regards to the capacitor 
Ci . For the results shown in Table U Ci was fixed to 
its known value while other parameters were calculated. 
In principle, Ci could also be obtained by optimization, 
but when we attempted this it was found that for this 
particular problem, shifting C'l away from its true value 
has only a very weak influence on the optimal cost value. 
Furthermore, some deficiency in the model causes the op- 
timization to prefer a drastically incorrect value for Ci 
while maintaining a correct value for the series equiva- 
lent capacitance, Ceq = 1/(1/Ci -I- I/C2). It seems that 
Ceq has matters more than Ci or C2 separately, which 
is reasonable because the fundamental frequency is given 
by /o = 1/(2^710;). 

To see if we could reduce or eliminate these effects, we 
tried adding some additional parameters to our model 
equations using a more accurate model of the transis- 
tor known as the Gummel-Poon model [i^. The full 
model involves an additional 10 parameters (excluding 



capacitive effects which we assume are negligible in our 
application). These include modeling internal resistance 
between the emitter of the transistor on the chip and the 
package lead which we connect to the rest of the circuit. 
There is a similar resistance between the collector and its 
package lead. The base resistance is assumed to vary with 
base current according to a particular formula including 
three adjustable parameters. There is also an assumed 
roUoff to the current gain involving one parameter. The 
Early voltage parameter compensates for changes in base 
thickness with collector voltage. 

It turns out that a significant improvement in the cost 
function can be achieved by merely adding the emitter 
resistance into the model equations. Because of the pres- 
ence of this resistor, —Ve is no longer the true base- 
emitter voltage, but instead it has a slightly higher volt- 
age —Ve- The equations for the transistor currents 
(Eqs. (g]) and (O) are now functions of Ve- The modi- 
fied voltage Ve is a function of Ve which must be found 
numerically by finding the solution to: 

VE^REls{l + ^)exp(^) +Ve- (21) 
Pf \ J 

Using only Re in addition to our previous parameters, 
the rms error in matching the data was reduced from 
about 0.0146 to about 0.0102 in volts. In addition, 
the value of /3f increased from 72 to 179, i.e. remark- 
ably close to the measured value of 180. We also tried 
other parameters from the Gummel-Poon model in var- 
ious combinations including the full set, but found no 
significant additional improvement in the cost function. 
Although the value of Pe was improved, the problem 
with Ci remains, so perhaps the Gummel-Poon model 
still retains significant deficiencies. These results were 
obtained using our first method (initial value), and are 
summarized in Table IIIII The parameters where also 
used to make a prediction of the fixed point, shown in 
Table ini 

D. Comparison of Models 

In order to better understand the effect of coupling 
on the optimization process, we again looked at how the 
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FIG. 13; Optimized rms error as a function of fixed, time inde- 
pendent coupling strength, (a) Real data analyzed with stan- 
dard model, (b) Real data analyzed with improved model, 
(c) Simulated data generated from standard model plus noise, 
and analyzed with standard model, (d) Simulated data gener- 
ated from improved model plus noise, and analyzed with im- 
proved model, (e) Simulated data generated from improved 
model, and analyzed with standard model. 



coupling strength effects the value of the rms synchro- 
nization error in fitting the data. These are shown in 
Fig. [131 Curve (a) (taken from Fig. [6] and included here 
for reference) shows the case where we are using exper- 
imental data and use our standard model (Eqs. (1-5)), 
while (b) shows this data represented using our improved 
model that include the effects of an internal emitter resis- 
tance Re- Note that in both cases the seemingly logical 
assumption that increasing the coupling strength will al- 
ways decrease the rms error is violated. Note that curve 
(b) starts at a significantly lower error than (a), a strong 
indication that it is a better model. In curves (c) and 
(d) we are looking to see whether the anomalous behav- 
ior of (a) and (b) are due to inaccuracies in the model 
or whether the presence of additive noise can produce 
a similar effect. Here we used simulated data with the 
same parameter values found in the experiment to which 
Gaussian noise was added. In case (c) the simulation 
used Eqs. (1-5), i.e. it did not involve the use of an 
internal emitter resistor, while in (d) that resistor was 
included. In both cases the noise level was adjusted to 
produce about the same rms error for m = as is found 
in the curves (a) and (b). Note that presence of noise 
can produce a positive slope effect, but seems to do so 
in a much less dramatic fashion than curve (a). Unlike 
curves (a) and (c), curves (b) and (d) seem much more 
similar to each other. This is perhaps an indication that 
the model used in those cases is a much better one. In 
the final curve (e), the simulated data was generated us- 
ing Re, but then it was analyzed with the model that 
lacks this resistor. No noise was added in that case, so 
the rms error is entirely due to a lack of ability of the 
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FIG. 14: Optimized coupling strength u{t) as a function of 
time with rj = 1.0. (a) Jagged curve from experimental data, 
analyzed without emitter resistance, (b) Smooth curve for 
simulated data with parameters matching the experimental 
data that was generated using an emitter resistance but ana- 
lyzed without (imperfect model). 



model to fit the simulated data. Note the similarity to 
curve (a), both peak at u approximately 1.7. It starts 
at a lower rms error than (a) because the real data does 
include some added noise. 

Model imperfections are yet another cause of the struc- 
ture seen in the time dependent coupling function u(t). 
We can demonstrate this by reanalyzing the simulated 
data of case (e) above, using our constrained method 
with time dependent coupling. The result shows much 
smoother features in u(t) than are generated by noise in 
the data. In Fig. [T3] we show just a portion of the u{t) 
function showing an overlay of the result for the simu- 
lated data with that of the experimental data. Note that 
the experimental data appears to show the same feature 
but with noise superimposed. This result suggests an- 
other way to identify imperfection in a model is 
to look for these larger scale smooth features in 
uit). 

IV. DISCUSSION 

We have developed a method, based on synchroniza- 
tion, for finding unknown parameters and predicting the 
future state of an actual chaotic experimental system, 
from limited measurements. This problem is difficult for 
several reasons: the dynamics are chaotic, there are many 
parameters and state variables to search over, we have 
only very limited measurements of the experimental sys- 
tem, and most importantly, we do not have an exact 
description of the experiment in our models. 

Because of the chaotic nature of the dynamics, the cost 
function is an especially rough surface in the space of 
parameters and initial state values of the response sys- 
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tern. This makes it difficult to locate the global mini- 
mum, which represents the point where the model sys- 
tem parameters match the data system parameters. This 
problem can be made dramatically easier by adding a 
term to the dynamics of the model system which couples 
the response system to the driver system. For sufficient 
coupling strength, this often has the effect of changing 
the largest conditional Lyapunov exponent from positive 
to negative, and stabilizes the synchronization manifold. 
This effect can be observed by noticing that the coupling 
smoothes out the cost function (Fig. S]). 

Even with the coupling, it is still difficult to search 
for a global minimum in a space with many dimen- 
sions. To perform the search we used a few different 
well-established optimization methods [s^, [H, [39| . This 
is mainly just a practical problem of efficiency. We did 
this with time independent coupling and time dependent 
coupling with almost identical results. Not surprisingly, 
the process takes much longer with time dependent cou- 
pling, because there are many more degrees of freedom 
to be optimized. 

We explored two different types of coupling. The sim- 
pler type is time independent coupling, which works quite 
well on our test problem. If the time series being ana- 
lyzed is short enough then coupling can be reduced all 
the way to zero at the end of the process. Otherwise 
it is reduced to the minimum level needed to maintain 
synchronization. We found an interesting effect where 
the optimal rms error in fitting the model to the data 
initially increases as the coupling strength is increased. 
This was very unexpected since the coupling is pulling 
the model variable towards the data. As the coupling 
is further increased, the error reaches a peak and then 
begins to decrease, so that for large coupling the behav- 
ior is what would be expected. The result is likely an 
artifact of the noise in the data, i.e. the coupling term 
is sometimes responding more to this noise than to the 
actual synchronization error. 

The other type of coupling is with time dependent u(t) 
which is optimized like the other parameters to yield the 
lowest cost function. This makes the problem more com- 
plicated, but gives us some additional, possibly useful, 
information such as the peaks in the residual u{t) which 
are robust and are localized in phase space. This struc- 
ture was also shown to be highly dependent on how u(t) 
is included in the cost function - a term that is quadratic 
in u leads to a very rich structure that mirrors the sen- 
sitivity of the synchronization error to the coupling as 
a function of time; a term that is linear in u leads to a 
few very sharp narrow peaks in the residual u{t), with 
the vast majority of time-steps having a u that is exactly 
zero. We suggest that this form of coupling may be best 
for maintaining synchronization for long data sets while 
causing minimal effect on the dynamics. 

By doing tests with simulated data we saw that there 
are three separate causes that contribute to the structure 
of the residual coupling u{t). 

• One is due to inaccuracies in the model. In re- 



gions where the model is inaccurate, the coupling 
increases to pull the model toward the data. This 
information could be useful for finding out where 
the model needs improvement. 

• The second is due to noise. We saw that u{t) 
changes significantly if only the particular sequence 
of random numbers used to generate the noise is 
changed, but there are persistent features. 

• The third cause is related to the direction of the pri- 
mary Lyapunov vector relative to the direction as- 
sociated with the coupling variable. This becomes 
relevant, even with a perfect model, when the time 
series is long enough to require coupling to coun- 
teract positive Lyapunov exponents. 

We also described a new method of coupling that is 
particularly useful for sparse data, i.e. data that has a 
time-step that is much too large to be used as the time- 
step for the model because it would lead to significant nu- 
merical error. This coupling method is applied discretely 
at the exact locations in time where the data points re- 
side. 

We considered the important effect of imperfect mod- 
els. In the Colpitts case, most of the circuit components 
behave in a well known nearly ideal way, and can be very 
accurately described by linear differential equations. The 
exception is the transistor. Even though it is moderately 
well described by the Ebers-MoU model, the slight in- 
accuracy that is present is sufficient to cause significant 
errors in the estimated values of some parameters. This 
can be verified by including an additional term in the 
equations - one that is part of the more detailed Gummel- 
Poon model [1^ . In many other cases that one might 
choose to study, the models are not as well established 
or as accurate and there will likely be similar problems 
but more severe. Our main interest in this simple case is 
that it provides a way of testing methods, that then can 
be applied to the harder cases. 

In this case we were then able to make improvements to 
the transistor model by adding the emitter resistance Re, 
but if we had no information besides Ve (t) we might have 
declared the original model good enough and stopped 
there. With that issue in mind, we then explored how 
the rms error of the fit depends on coupling strength 
(Fig. [HI). 

We proposed a test that can be done with just a mea- 
surement of one of the dynamical variables {Ve {tn ) in this 
case), and one prospective model. The main idea is that if 
the experiment is accurately described by the prospective 
model, then data simulated using the prospective model 
should be indistinguishable from the real experimental 
data. Here we are not comparing the real data directly 
to the simulated data, but instead comparing the results 
of the fitting process (performed using either the same 
model or a different model) applied to both real data and 
simulated data. We showed in Fig. [13] that when simu- 
lated data using the improved model (with noise added) 
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is fed into the fitting algorithm, the output (curve d) 
looks almost the same as the output of feeding the real 
data into the same fitting algorithm (curve b). Although 
there are many details of this test we do not understand 
yet, it seems to be useful way of checking the quality of 
a prospective model. 

Improvements to a model, when available, provide an- 
other test of our analysis - they are important when they 
cause a significant change in one or more of the estimated 
parameters or in the values of the variables that are in- 
accessible in the experiment. 

In a perspective on the ideas in this paper, we note 
that we addressed three problems: 

• when comparing data from an experimental source 
and a model where the time series has chaotic oscil- 
lations, one must regularize the search mechanism 
(minimizing a cost function) to avoid the irregular- 
ities in parameter and initial condition space so the 
impediments to the search of multiple local minima 
are removed, 

• one must select numerical methods that give the 
parameters and state variables over the observation 
period. 



• and one must develop good models to compare to 
the data, and be aware that seemingly small defects 
in a model can have a very strong impact on the 
results obtained. 

The methods we have described, including the two possi- 
ble choices for the second item that we explored, impact 
the quality of the model in an interesting fashion: the 
techniques can be used to explore the quality of a model, 
as we demonstrated in the investigation of two somewhat 
different models for the nonlinear transistor element in 
the Colpitts circuit. 
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